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In quantum chromodynamics (QCD) at nonzero chemical potential, the eigenvalues of the Dirac 
operator are scattered in the complex plane. Can the fluctuation properties of the Dirac spectrum 
be described by universal predictions of non-Hermitian random matrix theory? We introduce an 
unfolding procedure for complex eigenvalues and apply it to data from lattice QCD at finite chemical 
potential fi to construct the nearest-neighbor spacing distribution of adjacent eigenvalues in the 
complex plane. For intermediate values of n, we find agreement with predictions of the Ginibre 
ensemble of random matrix theory, both in the confinement and in the deconfinement phase. 

PACS numbers: 5.45.Pq, 12.38.Gc 



Physical systems which are described by non-Hcrmi- 
tian operators have recently attracted a lot of attention. 
They are of interest, e.g., in dissipative quantum chaos 
y, disordered systems ||,^, neural networks Q, and 
quantum chromodynamics (QCD) at finite chemical po- 
tential Q. Often they display unusual and unexpected 
behavior, such as a delocalization transition in one di- 
mension P]. Consequently, analytical efforts have been 
made to develop mathematical methods to deal with non- 
Hermitian matrices, see, e.g., Refs. [p^pOl. 

Of particular interest in the analysis of complicated 
quantum systems are the properties of the eigenvalues 
of the Hamilton operator. In the Hermitian case, it has 
been shown for many different systems that the spectral 
fluctuations on the scale of the mean level spacing are 
given by universal predictions of random matrix theory 
(RMT) pi. In QCD, one studies the eigenvalues of the 
Dirac operator, and it was demonstrated in lattice simu- 
lations that at zero chemical potential, the local spectral 
fluctuation properties in the bulk of the spectrum are 
reproduced by Hermitian RMT both in the confinement 
and in the deconfinement phase |12|] . 

If one considers QCD at nonzero chemical potential, 
the Dirac operator loses its Hermiticity properties so that 
its eigenvalues become complex. The aim of the present 
paper is to investigate whether non-Hermitian RMT is 
able to describe the fluctuation properties of the complex 
eigenvalues of the QCD Dirac operator. The eigenvalues 
are generated on the lattice for various values of fi. We 
define a two-dimensional unfolding procedure to sepa- 
rate the average eigenvalue density from the fluctuations 
and construct the nearest-neighbor spacing distribution 
P{s) of adjacent eigenvalues in the complex plane. The 
data are then compared to analytical predictions of non- 
Hermitian RMT. 

We start with a few definitions. At /i 7^ 0, the QCD 
Dirac operator on the lattice in the staggered formulation 
is given by [0| 
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with the link variables U and the staggered phases rj. 
The lattice constant is denoted by a, and color indices 
have been suppressed. 

We consider gauge group SU(3) which corresponds to 
the symmetry class of the chiral unitary ensemble of 
RMT |1J] . At zero chemical potential, all Dirac eigenval- 
ues are purely imaginary, and the nearest-neighbor spac- 
ing distribution P{s) of the lattice data agrees with the 
Wigner surmise of Hermitian RMT, 



Pw{s) 



32 



s2e-4«/^ 



(2) 



both in the confinement and in the deconfinement phase 
|l2| . This finding implies strong correlations of the eigen- 
values. In contrast, for uncorrclated eigenvalues P{s) is 
given by the Poisson distribution, Pp(s) — e^*. 

For fj. > 0, the eigenvalues of the matrix in Eq. (|l|) 
move into the complex plane. If the real and imaginary 
parts of the strongly correlated eigenvalues have approx- 
imately the same average magnitude, the system should 
be described by the Ginibre ensemble of non-Hermitian 
RMT 1^. As is the case in Hermitian RMT, we assume 
that the chiral structure of the problem does not affect 
the spectral fluctuation properties in the spectrum bulk. 

For a complex spectrum, we define P{s) to represent 
the spacing distribution of nearest neighbors in the com- 
plex plane, i.e., for each eigenvalue zq one identifies the 
eigenvalue zi for which s — \zi — Zf)\ is a minimum [n|. 
Clearly, this prescription is not unique, but it is the most 
natural choice and, more importantly, the only choice 
for which analytical results are available. After ensem- 
ble averaging, one obtains a function P{s, zq) which, in 
general, depends on zq. However, the dependence on zq 



can be eliminated by unfolding the spectrum, i.e., by ap- 
plying a local rescaling of the energy scale so that the 
average spectral density is constant in a bounded region 
in the complex plane and zero outside. This will be dis- 
cussed below. After unfolding, a spectral average over 
zq yields P{s). Of course, the question of whether, for a 
given data set, the average behavior of the spectral den- 
sity can indeed be separated from the fluctuations has to 
be answered empirically. 

In the Ginibre ensemble, the average spectral density 
is already constant inside a circle and zero outside, re- 
spectively H. In this case, unfolding is not necessary, 
and P{s) is given by [Q 
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where e„(a;) = X^m^o-^™/"^' ^^'-^ '^ ~ lo dssp{s) — 
1.1429.... This result holds for strongly non-Hcrmitian 
matrices, i.e., for Re(z) w Im(z) on average. In the 
regime of weak non-Hermiticity 0], where the typical 
magnitude of the imaginary parts of the eigenvalues is 
equal to the mean spacing of the real parts, the RMT 
prediction deviates from Eq. (|^). (Our definitions differ 
from those of Ref. |7| by a factor of i so that in our case 
it would be more appropriate to speak of weak non-anti- 
Hermiticity.) We shall comment on this regime below. 
For uncorrelated eigenvalues in the complex plane, the 
Poisson distribution becomes M 
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We now introduce a simple method to unfold complex 
spectra. This method will then be applied to the Dirac 
spectrum on the lattice at various values of /i 7^ 0, and 
the resulting P{s) will be compared to Eqs. (||) and (||). 
Assuming that the spectral density has an average and a 
fluctuating part, p{x,y) — pav(a;,y) +/3fl(2;,y), we need 
to find a map 



z' — x' + iy' = u{x, y) + iv[x, y) 



(6) 



such that pav(a;', J/') = 1. Conservation of the probability 
implies that p^vix' ,y')dx'dy' — dx'dy' — pii^{x,y)dxdy. 
Hence, pa.v{x,y) is the Jacobian of the transformation 
from {x,y) to {x',y'), 



Pa.vix,y) = \dxudyV - dyudxv\ 



(7) 



The left-hand side of Eq. (M) is given by the data. There 
are infinitely many ways to choose the functions u and 
V in Eq. (0) so that Eq. (M) is satisfied. (In general, 
a conformal map m) which fulfills Eq. (R) can only be 



constructed if p^vix^y) satisfies special conditions.) We 
choose y' = v{x,y) = y which yields pa,-v{x,y) = \dxu\ 
and, thus. 



x' = u{x,y) 



dtpa^{t,y) . 



(8) 



This corresponds to a one-dimensional unfolding in strips 
parallel to the real axis. For a fixed bin in y, p^vix, y) is 
obtained by fitting p{x, y) to a low-order polynomial. 

Let us now discuss how the lattice data were obtained. 
The simulations were done with gauge group SU(3) on a 
6'^ X 4 lattice using fi = 6/(7^ = 5.2 in the confinement re- 
gion and /3 — 5.4 in the deconfinement region for Nf = 3 
flavors of staggered fermions of mass ma = 0.1. For 
each parameter set, we sampled 50 independent configu- 
rations. The gauge field configurations were generated at 
/U = 0, and the chemical potential was added to the Dirac 
matrix afterwards. This procedure requires an explana- 
tion, since it is known that the quenched approximation 
at /^ 7^ is unphysical g. However, despite serious ef- 
forts ||l^ there is currently no feasible solution to the 
problem of a complex weight function in lattice simula- 
tions. In a random matrix model, the statistical effort 
to generate configurations including the complex Dirac 
determinant was shown to grow exponentially with p^N , 
where N is the lattice size [ [l6| . 

Typical eigenvalue spectra are shown in Fig. n^ for 
four different values of p (in units of 1/a) at /3 = 5.2. 
As expected, the size of the real parts of the eigen- 
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FIG. 1. Scatter plot of the eigenvalues of the Dirac oper- 
ator (in units of 1/a) in the complex plane at various values 
of fj, for a typical configuration (generated at /x = 0) in the 
confinement region at /3 — 5.2. 



values grows with fi, consistent with Ref. |lj]. Since 
the average spectral density is not constant, we have 
to apply the unfolding method defined above. Figure g 
shows the effect of the unfolding procedure on the eigen- 
value density in the complex plane. P{s) is then con- 
structed from the unfolded density and normalized such 
that /p°° ds s P{s) = 1. 

We have performed several checks of our unfolding 
method, (i) If only parts of the spectral support are con- 
sidered, the results for P{s) do not change. This means 
that spectral ergodicity holds, (ii) If the spectral density 
has "holes" (see Fig. |l| for fi — 1.0 and 2.2), we split the 
spectral support into several pieces and unfold them sep- 
arately. This is justified by spectral ergodicity. (iii) Un- 
folding each spectrum separately and ensemble unfolding 
yield identical results for P(s). (iv) The results for P{s) 
are stable under variations of the degree of the fit poly- 
nomial and of the bin sizes in x and y. Thus, we can 
conclude that it was indeed possible with our unfolding 
method for complex spectra to separate the average spec- 
tral density from the fluctuations, which is the prerequi- 
site for comparisons with analytical RMT predictions. 

Our results for P{s) are shown in Fig. 0. There are 
minor quantitative but no qualitative differences between 
confinement and deconfinement phases, which is consis- 
tent with our findings at /x = (second reference of 
|p2[). As a function of /i, we expect to find a transition 
from Wigner to Ginibre behavior in P{s). This is indeed 
seen in the figures. For fj, = 0.1, the data are still very 
close to the Wigner distribution of Eq. (0) whereas for 
0.5 < /i < 0.7 {n = 0.7 not shown) we observe nice agree- 
ment with the Ginibre distribution of Eq. (0) . Values of 
fi in the crossover region between Wigner and Ginibre be- 
havior (0.1 < fi < 0.3) correspond to the regime of weak 
non-Hermiticity mentioned above. In this regime, the 
derivation of the spacing distribution is a very difficult 
problem, and the only known analytical result is P{s, zq) 
for small s, where zq is the location in the complex plane 
(i.e., no unfolding is performed) 0. The small-s be- 
havior of Eqs. (g) and (g) is given by Av(s) 
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FIG. 2. Eigenvalue density of the Dirac matrix in the com- 
plex plane for a typical configuration at /3 = 5.2 and /i = 0.7. 
The left plot shows the eigenvalues before unfolding, the right 
plot after unfolding, respectively. 
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FIG. 3. Nearest-neighbor spacing distribution of the Dirac 
operator eigenvalues in the complex plane for various values 
of fj, in the confinement (/3 — 5.2, left) and deconfinement 
{13 = 5.4, right) phase. The histograms represent the lattice 
data. The solid curve is the Ginibre distribution of Eq. (H), 
the short-dashed curve in the upper two rows the Wigner 
distribution of Eq. (H), and the dotted curve in the lowest 
row the Poisson distribution of Eq. (H), respectively. 



Pg(s) oc s^, respectively, and in the regime of weak non- 
Hermiticity we have P{s, zq) oc s" (for s ^ 1) with 
2 < a < 3 0. This smooth crossover from a = 2 to 
a = 3 is also observed in our unfolded data. 

For fj. > 0.7 the lattice results for P{s) deviate sub- 
stantially from the Ginibre distribution. It is also ev- 
ident from the appearance of the spectra for /x = 1.0 
and 2.2 in Fig. fi^ that the global spectral density of the 
lattice data is very different from that of the Ginibre en- 
semble. While this docs not immediately imply that the 
local spectral fluctuations are also different, it is an indi- 
cation for qualitative changes. In particular, the results 
for fj, = 2.2 in Fig. y could be interpreted as Poisson 
behavior, corresponding to uncorrelated eigenvalues. In 
the Hermitian case at finite temperature, T, lattice sim- 
ulations show a transition to Poisson behavior only for 
/3 — > oo when the physical box size shrinks and the theory 
becomes free |12|] . Here, however, we have to remember 
that the chemical potential was neglected in the gener- 
ation of the equilibrium gauge fields so that conclusions 
about the physical content of this observation should be 
drawn with care. Specifically, we are not entitled to make 
definite statements on possible connections between the 
finite-density phase transition (expected at smaller val- 
ues of fi) and the deviations from Ginibre behavior. A 
more plausible explanation of the transition to Poisson 
behavior is provided by the following two (related) ob- 
servations. First, for large /i the terms containing e^ in 
Eq. (Wj dominate the Dirac matrix, giving rise to uncor- 
related eigenvalues. Second, for // > 1.0 the fermion den- 
sity on the 6^ X 4 lattice reaches its maximum value given 
by the Pauli exclusion principle. In any event, the mech- 
anisms behind this transition are interesting and deserve 
further study. 

If a solution to the problem of the generation of config- 
urations with a complex weight function becomes avail- 
able, it would be very interesting to redo our analysis. 
Since previous computations in the deconfinement phase 
at /i = have verified the RMT predictions [Q, and 
since the present simulations at /3 = 5.4 are already in 
the deconfinement phase for all values of fi, we expect 
the observed Ginibre behavior to persist. Of course, this 
should be verified by a full finite-density simulation. 

In conclusion, we have proposed a general unfolding 
procedure for the spectra of non-Hermitian operators. 
This procedure was applied to the QCD lattice Dirac 
operator at finite chemical potential. Agreement of the 
nearest-neighbor spacing distribution with predictions of 
the Ginibre ensemble of non- Hermitian RMT was found 
between /i = 0.5 and /i = 0.7 in both confinement and 
deconfinement phases. The deviations from Ginibre be- 
havior for smaller values of /i are well understood whereas 
the deviations for larger values of fi toward a Poisson 
distribution require a better understanding of QCD at 
finite density. In this context, an interesting observation 
is that the results for P{s) in the non-Hermitian case are 



rather sensitive to fj, whereas they are very stable under 
variations of T in the Hermitian case ||lj] . It would be 
interesting to apply our unfolding method to other non- 
Hermitian systems [Q,Q . 
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